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ABSTRACT 

Galaxies selected at 170 by the ISO FIRBACK survey represent the brightest ^ 10% of 
the Cosmic Infrared Background. Examining their nature in detail is therefore crucial for con- 
straining models of galaxy evolution. Here we combine Spitzer archival data with previous 
near-IR, far-IR, and sub-mm observations of a representative sample of 22 FIRBACK galaxies 
spanning three orders of magnitude in infrared luminosity. We fit a flexible, multi-component, 
empirical SED model of star-forming galaxies designed to model the entire ^ 1 - 1000 fim 
wavelength range. The fits are performed with a Markov Chain Monte Carlo (MCMC) ap- 
proach, allowing for meaningful uncertainties to be derived. This approach also highlights 
degeneracies such as between and /3, which we discuss in detail. From these fits and stan- 
dard relations we derive: iiR, Lpah, SFR, ry. A/*, Mdust, 7d, and [3. We look at a variety 
of correlations between these and combinations thereof in order to examine the physical na- 
ture of these galaxies. Our conclusions are supplemented by morphological examination of 
the sources, and comparison with local samples. We find the bulk of our sample to be consis- 
tent with fairly standard size and mass disk galaxies with somewhat enhanced star-formation 
relative to local spirals, but likely not bona fide starbursts. A few higher-z LIGs and ULIGs 
are also present, but contrary to expectation, they are weak mid-IR emitters and overall are 
consistent with star-formation over an extended cold region rather than concentrated in the 
nuclear regions. We discuss the implications of this study for understanding populations de- 
tected at other wavelengths, such as the bright 850 ftm SCUBA sources or the faint Spitzer 
24 fim sources. 

Key words: galaxies: fundamental parameters, galaxies: starburst, infrared: galaxies, submil- 
limetre 



1 INTRODUCTION 

Understanding the full infrared spectral energy distribution (SED) 
of galaxies is essential for a complete picture of star-formation in 
the Universe. Measurements of the cosmic background radiation 
allow us to infer that about half the energy ever g enerated by stars 
was r eprocessed by dust into the infrared (see iHauser &Dweld 
I2OOII for a review). This emission is increasingly important at 
higher redshifts where the star-formation density of the Universe 
is larger than today. Modelling the contribution of different galaxy 
populations to the Cosmic Infrared Background (CIB) requires 
detailed knowledge of the SEDs of star-forming galaxies. The 
variation in SED shapes is a key uncertainty when comparing 
populations selected at different wavelengths or testing galaxy 
evolution models. 

In this paper, we discuss the full infrared SED (~ 1 - 1000 /im), 
which roughly spans the wavelengths between stellar and 
synchrotron-dominated emission. Traditionally, studying this 



entire range at once has been difficult, since mid-IR (<60/im) 
observations could not reach much beyond the local Universe^, 
while sub-mm obser vations also only ex ist for very local, IR- 
bright, galaxies (e.g. Dunne & Eales '' 200lll . or else the blank-sky 
Sub-mm Common U ser Bolometer Array (SCUBA) sources 
which peak at 2 ~ 2-3 iChapman et aljEoOSh . The latter typically 
have only one or two other wavelength detections in addition 
to the SCUBA (850 fim) one, which makes their in terpretation 
particularly dependent on the SED model assumed (see lBlain et alJ 
1200^ for a review). Due to these past observational limitations, 
we still do not know (beyond some generalized trends) the full 
range of galaxy SED shapes, how exactly they are related to the 
underlying physical conditions in the galaxy, and therefore how 
they may vary across cosmic time as the galaxies evolve. 

^ Except for the deepest ISOCAM 15 observations which peak at 
2 ~ 0.7 ( Rodighiero et al. 2004) 
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2 DATA 

The galaxies we focus on come from a sub-sample of FIRBACK 
sources with radio detections which were observed with SCUBA. 
The full sample of 30 targets was first presented in S03 where de- 
tails of the sample selection as well as the near-IR {K) and sub- 
mm observations and data reduction are given. In particular, in 803 
we addressed the question of whether or not the radio detection 
requirement introduces a bias and found that apart from potential 
HyLIGs (I/iR, > lO'^^ Lq) at z ~ 1.5, and spurrious sources, the ra- 
dio detection does not bias us with respect to the FIRBACK cat- 
alogue as a whole. Here we additionally use Spitzer archival ob- 
servations of the ELAIS-Nl field obtained as part of the Spitzer 
Wide-area IR Extragalactic Survey (SWIRE, Lonsdale et al. 2003). 
We focus on sources with unique identifications in the SWIRE data 
(see Section 2.2), which leaves us with a sample of 22 sources (see 
Table0. 



Figure 1. IRAC 3.6 /xm images for bright galaxies in our sample. Note 
the prevalence of disk-like morphologies with bright nuclei (Nl-004 is a 
spectroscopically confirmed weak AGN). Boxes here are approximately 
1.4 arcmin wide. For comparison, ISO's 170 fim beam is ~ 90 ai'csec. 



With the advent of the Spitzer Space Telescope iWerner et alJ 
l2004h we can for the first time observe the mid-IR properties 
of large numbers of sources over a cosmologically significant 
extent in redshift (e.g. Lonsdale et al. 2003). Spitzer covers the 
range 3-160 /im. The obvious next step toward characterizing 
the full infrared SEDs of galaxies is therefore to link Spitzer 
observations with longer wavelength samples, especially including 
sub-mm observations. In this regard, the quality and quantity of 
the available sub-mm data are the limiting factor. 
The sample discussed here is selected from the 1 70 um FIR- 
BACK (Far-IR BACKground) ELAIS-Nl catalog iPueetetal. 
Ll999; Dole et al. 2001). The selection is based on an existing 
radio detection, which we foUowed-up with both deep near-IR 
and sub-mm observations (Scott et al. 2000; Sajina et al. 2003, 
hereafter S03). We find that in terms of the mid-IR properties, the 
sample selection is largely unbiased with respect to the FIRBACK 
population as a whole. Our previous studies suggest that the 
sample consists primarily of 2 < 0.3 ordinary, spiral-like galaxies, 
rich in cold (r<40K) dust, with ~l/6 of the sample consist- 
ing of Ultraluminous Infrared Galaxies (ULIGs) at 2 0.5-1. 
Spectroscopic follo w-up of this and relat ed sub-samples support 
these findings (.Chapman et aljE'002l : IPatris et al...2003: Dennefeld 
I2OO5I hereafter D05). The higher-z fraction therefore represents a 
bridge population between the local Universe and distant, dusty 
star-formers such as the SCUBA blank-sky sources. In general, 
the importance of FIRBACK galaxies is that they represent the 
brightest galaxies at 170 /xm, contributing ~ 10% to the CIB and 
selected at a wavelength near the peak of the CIB spectrum. 
We use archival Spitzer observations of the ELAIS-Nl field in 
order to extend the known SEDs of the above sample into the 
mid-IR wavelength range. We fit these SEDs with a phenomeno- 
logical model motivated by different physical origins for the 
emission. These fits allow us to discuss the physical characteristics 
of our galaxies as well as trends within the sample. Details of the 
fitting procedure and some related issues are included in a set of 
appendices. 

Throughout the paper we assume a flat Universe with 
i7o =75 km s"^ Mpc"\ = 0.3, and Qa=OJ. 



2.1 Treatment of Spitzer data 

The basic data reduction and calibration is already performed on 
the individual frames obtained from the Spitzer archive. We use the 
Starlink software package CCDPACK to rescale, align and coadd the 
observations into common mosaic images for each of the 4 IRAC 
bands and the MIPS 24 fim band. Furthermore, the Starlink PHO- 
TOM and GAIA packages are used to perform aperture photometry 
on the sources. To ensure that the sky-annuli chosen fairly repre- 
sent the sky/background we monitor the sky values obtained for 
all of these and modify the annular region for any galaxy whose 
sky value is more than 2 a above the average sky. The errors are 
then computed from the sky variance. We find unambiguous coun- 
terparts for nearly all radio/SCUBA sources in S03 (the exceptions 
being NI-032, and Nl-034). 

Since we started on this proj ect, the SWIRE tea m has released their 
own catalogues of the field ( Surace et all20o4) . These allow us to 
double check our IRAC and MIPS 24 fim photometry, as well as 
add a 70 ^lm point where available. For Nl-004, NI-007, and NI- 
012 there are no 24fim or 70 /im fluxes due to missing data. In 
addition to the SWIRE 70 ^m catalogue, we extract fluxes for three 
faint sources: NI-040, NI-064, and NI-077. In all cases, the aper- 
ture resulting in the maximum flux was used in order to ensure all 
emission is accounted for. The only remaining source is Nl-078 
which is near the edge of the image and thus a confident flux can- 
not be obtained. 

This is the only source without any data points between 24 ^m and 
170 ^m. For all of Nl-004, Nl-007, and Nl-012 we have ISO- 
CAM 15 and IRAS 60 ^m fluxes compensating for the miss- 
ing MIPS data. For the few cases which have both a 60 /im (see 
D05) and 70 ^m detection (Nl-OOl, Nl-002, and Nl-016), we find 
that the 60 fim fluxes are somewhat higher than the 70 ^m ones 
contrary to any reasonable SED (given the S7o/Sno colours); but 
the difference is within the 20% calibration uncertainty assigned 
to the 70 ^im flux. The SWIRE catalogue 70 /im flux for NI-OOI 
(198 mjy) was most discrepant originally; however, it is an ex- 
tended source and we find that a flux of 233 mJy is more accurate. 
This 20% difference is the most severe we expect due to aperture 
effects for this sample. We do not explicitly use the few available 
160 pm points, but within the uncertainties they are consistent with 
the ISO 170 fim points. 

Near-IR y-band data for about half of our sample are available from 
the band-merged ELAIS catalog teowan-Robinson et aljEo04l) . 
while the A^-band data come from our previous work (S03). To con- 
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Table 1. Multiwavelength data for our sample. All errors are 1 a estimates. 



Source* 


,91 ■3°' 


Q b 
'-'2.2 


'-'3.6 


'-'4.5 




^8.0 


o a 

024 


■^70 


■^170 


■^450 


^850 


^ 1 .4GHz 




mJy 


mJy 


mJy 


mJy 


mJy 


mJy 


mJy 


mJy 


mJy 


rajy 


mJy 


mJy 


N 1-001 


6 28 


1 ^1 -|-0 54 


2 97^0 37 


1 98±0 26 


4 67-|-o 4 


12 25iO 66 


19 97i 1 71 


233 


597i72 


-3 0± 14 


6 1 ± 1 6 


74±0 23 


N 1 -002 


4 58 


5 55zbO 41 


5 35zbO 43 


3 7g-|-0 37 


9 6±0 59 


23 97iO 92 


(J g7-j- ] 21 


1 12 


544 J; 69 


14 4± 12 4 


4 4± 1 1 


64±0 04 


N 1-004 


9 16 


J 3] -|-0 14 


5 38zbO 44 


3 7g-|-o 36 


6 05 ±0 47 


1 1 42it0 64 






391it58 


32 5±7.2 


3 6± 1 4 


88±0 13 


N 1 -007 


4 40 


3 5zbO 26 


3 65zb0 36 


2 41 J^o 29 


4 09±0 38 


1 1 87iO 65 






338it54 


23 4±8 1 


4 4± 1 6 


1 04±0 12 


N 1 -009 


9.63 


57_j_Q 2 


6.02zb0.46 


3.96±0.37 


7. 05 ±0.5 


14.03it0.7 


4.26±0.8 


96 


313±52 


10.6±7.6 


3.5±1.5 


1.15±0.1 1 


N 1-0 1 2 


I 75 


[ 84-|-0 3 


1 15zbO 2 


Q 75-|-o 16 


81±0 17 


4 24i0 39 






302it51 


9 2± 10 


1 5± 1 6 


31 ±0 07 


Nl-On 




Q n-|-o 01 


16zb0 07 


Q i4-|-o 07 


18±0 08 


44i0 12 


2 21 itO 58 


83 


294ib5 1 


18 8±9 9 


0± 1 5 


52±0 15 


N 1-0 1 5 


36 


80zbO 06 


Q 57-|-0 14 


53±0 14 


Q 37-|-o 1 1 


2 05ibO 27 


3 16it0 68 


49 


294i51 


-3 4±7 7 


1 4± 1 6 


52±0 07 


N 1-0 1 6 


2 71 


3 5zbO 26 


I 78-|-0 25 


y 41-|-{) 22 


1 61 ±0 24 


7 6it0 52 


6 81 i 1 


160 


289it50 


34 8± 16 7 


1 5± 1 2 


1 55±0.13 


N 1 -024 


1 14 


y 3g-|-0 03 


J i5-|-o 2 


Q 77-|-o 16 


Q g-|-Q [§ 


5 4^0 44 


4 38itO 81 


108 


266it49 


32 3±7 5 


2 9± 1 3 


75±0 02 


Nl-0''9 


I 10 


1 27zb0 09 


1 24zb0 21 






4 ^5itO 39 


4 3itO 8 


72 


229it46 


90 14 1 


5± 1 7 


n 6Q^n 05 


Nl-031 


1.93 


2.65±0.19 


1.45 ±0.22 


0.91±0.18 


1.44±0.23 


3.53±0.35 


541±0.9 


62 


225 ±46 


9.2±13.2 


1.9±1.1 


0.43±0.06 


N 1-039 




0.32±0.05 


0.28±0.1 


0.25±0.09 


0.19±0.08 


0.8±0.17 


1.3±044 


44 


205 ±44 


I0.9±86.8 


-0.1±2.3 


0.58±0.07 


N 1-040 




O.OlzbO.OI 


0.08±0.05 


0.07±0.05 


0.06±0.05 


0.05±0.04 


0.69±0.32 


26 


205 ±44 


29.2±20.5 


54±1.1 


0.33±0.03 


Nl-041 


0.38 


0.88±0.06 


0.66±0.15 


0.5±0.13 


0.58±0.14 


4.14±0.38 


3.89±0.76 


75 


204±44 


20.4± 156.7 


-0.1±2.5 


0.76±0.06 


N 1-045 


1.39 


1.16±0.02 


1.01±0.19 


0.66±0.15 


0.69±0.16 


3.29±0.34 


2.23±0.58 


64 


198±44 


15.3±8.3 


3.0±14 


0.43±0.06 


N 1-064 




0.03±0.0I 


0.15±0.07 


0. 11 ±0.06 


0.12±0.07 


0.09±0.06 


1.38±046 


14 


166±42 


35.2±13.9 


5.1±1.2 


0.23±0.04 


N 1-068 


0.31 


0.51 ±0.04 


0.3I±0.1 


0.31±0.11 


0.28±0.1 


1.98±0.26 


3.96±0.76 


62 


165 ±42 


15.1±7.6 


2.2±14 


0.44±0.05 


N 1-077 


0.31 


0.42±0.07 


0.37±0.11 


0.29±0.1 


0.23±0.09 


1.4±0.22 


1.29±044 


42 


159±41 


5.9±7.3 


l.l±1.3 


0.40±0.10 


N 1-078 




0.04±0.0I 


0.13±0.07 


O.I ±0.06 


0.09±0.06 


O.I3±0.07 


0.6±0.3 




158±41 


35.2±8.7 


5.7±1.3 


0.24±0.04 


N 1-083 




0.46±0.08 


0.3 1 ±0.1 


0.25±0.09 


O.I8±0.08 


0.68±0.16 


1.21 ±043 


43 


150±41 


16.2±16.0 


0.7±1.2 


0.55±0.03 


Nl-101 


0.38 


0.55±0.09 


0.58±0.14 


0.38±0.12 


0.46±0.13 


2.34±0.29 


2.66±0.63 


46 


136±40 


19.8±7.5 


0.9±1.5 


0.39±0.05 



* The naming scheme follows Dole et al. (2001). 

J-band magnitudes from Rowan-Robinson et al. (2004) where 20% uncertainty is assumed. 
UKIRT and SCUBA fluxes from Sajina et al. (2003). 
^ IRAC and MIPS fluxes from archival SWIRE observations, this work; S-jq from the SWIRE catalogue, all assumed to have an uncertainty of 20%. 
''■ ISOPHOT 170 ;jm data from Dole et al. (2001). 
' VLA 21 cm fluxes from Ciliegietal. (1999). 



vert to flux densities, we used zero points of 1600 Jy and 670 Jy for 
the /- and A'-bands respectively. 

The Spitzer fluxes are presented in TableQalong with the rest of the 
multiwavelength data used here. A few of the sources have some (or 
all of) U, G, R, ISOCAM 15 fim, IRAS 60 ^m and 100 /^m fluxes, 
which are given in DOS. In our present study, these are used primar- 
ily for consistency checks. 

2.2 Possibility of multiple sources 

We use the SWIRE ELAIS-Nl catalog to further investigate the 
nature of the FIRBACK sources and in particular the relation of 
our radio-selected sources to the full sample. The 90arcsec ISO 
beam means that the chance of multiple sources contributing to the 
170 ^m flux is not expected to be negligible. We perform a simple 
test of this effect by taking every 24 /im source within a 45 arcsec 
radius of the ISO position. We then add all these fluxes and com- 
pare with the single brightest source within the centroid. The re- 
sults are presented in Fig.|2| where we see that indeed nearly half 
of all sources regardless of whether the full or 4(t catalog 
is used are not well described by assigning a single counterpart. 
This is a useful qualitative test, although it is unlikely to be cor- 
rect in detail, because of complications such as foreground AGB 
stars, or combinations such as a bright local galaxy near a some- 
what higher-2 LIG/ULIG (where the latter might still contribute 
the bulk of the 170 ^m flux). Nevertheless, it is evident that nearly 
half of the 170 /.im-selected sources have 2 or 3 roughly equivalent 
mid-IR counterparts. Considering the position of the counterparts 
within the beam does not appear to affect this significantly. This 
fraction drops slightly for our radio-selected sample, although it is 
still significant. This general conclusion is in agreement with the 
results of Dennefeld et al. (2005) who find 28/56 sources in the 4 a 
catalog to have firm, unique identifications. Here we add that this 
fraction is unlikely to change significantly for the 3 a catalog. We 
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Figure 2. Quantitative study of the potential effect of multiple sources 
within the ISO beam at 170 fim. These histograms show the distribution 
of flux differences at 24 /im between using just the brightest 24 fim source 
within the ISO beam versus using all sources within a 45 arcsec eiTor cir- 
cle (i.e. 100 (1 - S24,max/J2 '5'24))- The solid line is for aU FIRBACK 
sources, while the dashed Une is for the sources in the > 4 cr catalog. The 
dotted line is for our sub-sample of 30 radio-selected sources. 

also find similar results for our FIRBACK sub-sample of 30 targets, 
with 28% having fainter 24 fim sources in the ISO beam contribut- 
ing > 30% of the flux of the brightest source (dotted line in Fig.|2}. 
Because of strong ambiguities, we remove 8 FIRBACK sources 
from our sub-sample. We leave 2 ambiguous sources (Nl-015 and 
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Table 2. Derived properties for our sample. Errors are marginalized 68% confidence limits. 



Source 


z" 


off 


d 


(3 


log A/j 




log A/* 


^3 — 1000 


1 T" 
log -'jpAH 




Quality^ 








kpc 


K 




IVIQ 




lVi0 






M Q/yr 




Nl 


001 


030 


4 41 


25 2 it 4 7 


1 9 ± 5 


7 10 ± 22 


95 ± 10 15 


9 44 ± 09 


10 06 ± 06 


9 07 ± 2 1 


1 72 ± 29 


2 


Nl 


002 


064 


9 55 


21 9 lb 6 1 


2 1 ± 6 


7 62 ± 24 


5 4 1 ± 5 60 


10 47 ± 16 


10 62 ± 06 


10 15 ± 20 


5 3 ± 1 1 


2 


Nl 


004 


064 


14 29 


23 lb 5 6 


2 1 ± 5 


7 56 ± 19 


3 31 ± 3 50 


10 43 ± 13 


10 60 ± 08 


9 71 ± 24 


5 7 ± 1 2 


] 


Nl 


007 


0.061 


6.10 


24.1 lb 7.1 


1.4 ± 0.7 


7.50 ± 0.23 


5.40 ± 5.03 


10.25 ±0.18 


10.42 ± 0.1 1 


9.71 ± 0.21 


3.6 ± 1.0 


2 


Nl 


009 


053 


5 56 


23 3 lb 6 6 


7 1 ± 6 


7 18 ± 37 


7 67 ± 7 37 


10 32 ± 13 


10 72 ± 08 


9 61 ±0 18 


2 02 ± 47 




Nl 


012 


066 


6 88 


27 4 lb 8 3 


1 7 ± 6 


7 35 ± 37 


2 27 ± 1 68 


9 65 ± 1 1 


10 44 ± 13 


9 26 ± 25 


4 2 ± 1 4 




Nl 


013 


((\ 461 


18 63 


32 5 it 1 1 6 


2 ± 8 


8 73 ± 34 


7 81 ± 4 12 


10 79 ± 16 


12 1 1 ± 1 1 


10 64 ± 23 


210 ± 55 




Nl 


015 


234 


12 37 


23 7 lb 5 1 


2 1 ± 5 


8 50 ± 37 


2 85 ± 1 63 


10 56 ± 16 


1 1 39 ± 12 


10 40 ± 25 


36 ± 13 




Nl 


016 


092 


6 10 


31 8 lb 8 


1 6 ± 6 


7 44 ± 35 


1 73 ± 8 02 


10 25 ± 15 


10 86 ± 08 


10 29 ± 26 


1 1 ± 1 9 




Nl 


024 


086 


5 44 


24 2 lb 7 9 


1 8 ± 7 


7 73 ± 23 


3 15 ± 2 7 1 


9 95 ± 16 


10 60 ± 09 


9 76 ± 20 


6 1 ± 1 6 




Nl 


029 


144 


9 75 


26 6 it 16 9 


2 ± 8 


7 88 ± 59 


3 50 ± 2 74 


10 40 ± 14 


10 93 ± 12 


10 09 ± 23 


12 2 ± 4.5 




Nl 


031 


0.063 


4.29 


21.9it 15.5 


2.4 ±0.8 


7.22 ±0.52 


2.08 ± 2.03 


9.81 ±0.12 


I0.16±0.13 


9.20 ± 0.25 


2.07 ± 0.69 




Nl 


039 


0.269 


14.21 


24.5 ±9.4 


2.5 ±0.6 


8.44 ±0.70 


7.21 ± 5.34 


10.53 ± 0.22 


ll.37±0.16 


10.21 ±0.24 


35 ± 13 




Nl 


040 


0.450 


26.34 


26.5 ±7.9 


1.8 ±0.7 


9.10±0.21 


15.61 ±7.36 


10.42 ±0.34 


ll.83±0.13 


9.85 ±0.41 


104 ± 33 




Nl 


041 


0.120 


7.85 


24.7 ± 12.2 


2.3 ±0.7 


7.82 ± 0.77 


4.71 ±2.51 


10.01 ±0.14 


10.77 ±0.14 


9.93 ±0.23 


8.2 ± 3.0 




Nl 


045 


(0.18) 


15.74 


28.1 ± 14.6 


1.1 ±0.8 


8.08 ±0.57 


2.55 ± 0.00 


10.53 ±0.16 


11.08±0.I3 


10.23 ±0.16 


17.4 ±6.0 




Nl 


064 


0.910 


17.54 


31.1 ±6.3 


1.8 ±0.7 


9.32 ±0.21 


29 ± 30 


11.47 ±0.25 


12.45 ±0.13 


11.73 ±0.39 


410± 130 




Nl 


068 


(0.16) 


9.69 


26.7 ± 20.3 


2.1 ±0.8 


7.83 ±0.75 


3.41 ± 1.99 


9.99 ±0.15 


10.93 ±0.14 


9.88 ±0.20 


13.0 ±4.4 




Nl 


077 


(0.20) 


11.60 


27.1 ± 14.5 


2.0 ±0.7 


7.91 ±0.65 


3.25 ± 2.27 


10.I7±0.I5 


10.92 ±0.15 


10.06 ±0.21 


12.0 ±4.6 




Nl 


078 


(0.91) 


16.33 


69.8 ± 15.0 


1.1 ±0.5 


9. 16 ±0.25 


27 ±30 


11.49 ±0.26 


12.90 ±0.22 


11.59 ±0.38 


1360 ±570 




Nl 


083 


(0.31) 


14.17 


30.8 ± 18.3 


2.0 ±0.7 


8.21 ±0.58 


4.81 ±4.19 


10.65 ± 0.23 


ll.40±0.16 


10.34±0.19 


35 ± 15 




Nl 


101 


0.060 


6.16 


24.6 ± 17.8 


2.1 ±0.8 


7.07 ± 0.75 


4.94 ± 3.33 


9.36 ±0.18 


9.89 ±0.14 


9.01 ±0.19 


1.07 ±0.44 





^ Spectrcscopic redshifLs from D05 and Chapman et al. (2002) (Nl-040 and Nl-064); photometric redshifts given in brackets (see Appendix B). 

The observed half-light radius at 8 ^m converted to physical distance. 
^ The integral under our full PAH template. 

As per Kennicutt (1998). For consistency Ijg — 1000 ^^^^ here. 

Quality flag: l=good; 2=x^g,^ > 2; 3=t-y poorly fit (see Section A.l). 



Nl-039) because they have Zspcc ~ 0.2-0.3 and this test is incon- 
clusive for higher-z sources. This leaves 22 sources for which the 
correct counterpart is reasonably secure, with other possible coun- 
terparts contributing a probable amount to the 170 /^m flux which 
is of order the flux uncertainties or less. Since the beamsize at the 
other wavelengths is so much better than at 170 /xm, the effects of 
multiple source contributors to the other fluxes are negligible. The 
effects of multiple contributors are already partially included in the 
170 ^m uncertainties, which contain confusion noise (see Lagache 
& Puget 2000). Therefore after removing the 8 ISO sources which 
have the most ambiguous identifications at 24 /im, we are confi- 
dent that the effects of multiple counterparts are not significant. We 
concentrate on this new sub-sample of 22 sources in the rest of this 
paper. 



3 SED MODEL 

3.1 Interpreting the observed SEDs 

The observed SED of a galaxy depends on the properties of the 
stellar populations (e.g. ages and metallicities), the dust model 
(e.g. composition, size distributions), and the relative geometry 
of the two. The full range of possible combinations of all these 
is such that, fully accounting for all possible processes is todate 
not possible. Radiative transfer models where stars form inside 
giant molecular clouds and gradually disperse into the diffuse 
medium exist iSil va et al. 1998; Efstathiou & Rowan-Robinson 
l2003t l Popescu et alj 2000), but tend either to have too many free 
parameters, or not cover the full range in SEDs observed. Self- 
consistent models of the stars and dust (i.e. energy a bsorbed equals 
energy emitted) have been done on local galaxies ( Galli ano et alj 
[2003); however, require exceptionally well sampled SEDs from the 
UV to the sub-millimeter. We choose a modelling approach (see 
Section 3.2) that is data-driven: i.e. our model is no more compli- 
cated than our data allow us to fit, yet is flexible enough to describe 
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Figure 3. An example of the SED model used here. This includes a grey- 
body {dotted line), a warm power-law (short-dash), PAH emission (long- 
dash), and unextincted stellar emission (dot-dash) with e^'^" extinction 
applied. The thick solid line is the total. 



the full range of SEDs observed in the sample. Such a simple phe- 
nomenological model has the effect of smoothing over the under- 
lying messiness of galaxy SEDs, while providing us with best-fit 
parameters (such as opacity and dust temperature) those values are 
some 'characteristic' or 'average' values of an underlying distri- 
bution. The approach is therefore well suited to broadband data, 
such as presented for our sample in the previous section, which is, 
in essence, spatially integrated over the entire galaxy, including its 
quiescent (i.e. cirrus) and potentially active (i.e. starburst) environ- 
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ments. 

The main limitation in interpreting our results in terms of the phys- 
ical nature of our galaxies is due to the fact that, as discussed in the 
previous section, the available optical/UV data (U, G, and R) for 
our sample are quite scarce and thus here we only discuss the spec- 
tral properties of our galaxies longward of and including the near- 
IR (~ 2 /im), which automatically precludes us from directly dis- 
cussing young stars. Thus our restriction of fitting only the ^ 2 jim 
spectra means that: 1) we cannot balance the energy absorbed in 
the optical/UV with the energy re-emitted in the infrared, 2) we 
cannot explicitly discuss either the star-formation history or age 
distribution of the stellar populations of our galaxies, and 3) we 
cannot constrain multiple optical depths (e.g. molecular cloud and 
cirrus). These limitations are largely addressed post-fitting, i.e. in 
the interpretation of the best-fit model parameters and how they 
are used to infer underlying physical characteristics. For example, 
we confine our stellar model to an old stellar population template 
(see Section 3.2), but when deriving stellar mass, we use mass-to- 
light ratios which include the effect of young stars on the /f-band 
(Section 4.3). We know that young stars form inside dense molecu- 
lar clouds and later disperse into the diffuse 'disk' environment. 
Thus the characteristic optical depth absorbing the young stars' 
power is likely not the optical depth we derive in the near-IR which 
is dominated by the older stellar populations. We address this in 
Section 5.2, where the dust masses are derived both from emis- 
sion (predominantly young stars) and absorption (predominantly 
old stars in our case). 

With the above philosophy in mind, we proceed with the model 
procedure outlined below. 

3.2 Details of Model Procedure 

As discussed above, our modelling approach is driven by the de- 
sire to have a phenomenological model which fits known nearby 
galaxy SEDs robustly and yet has the flexibility (i.e. number of 
phenomenological parameters only slightly less than the number of 
data points we are trying to fit) which allows us to fit a wide range 
of SED types, and in addition we want to be able to extract phys- 
ical parameters with meaningful uncertainties (in particular such 
that degeneracies between parameters can be highlighted). With 
these goals clearly in mind, our choice of detailed SED model is de- 
pendent on the quality and wavelength coverage of the data, while 
the statistical approach we use is the Markov Chain Monte Carlo 
method. 

We model the SED as a sum of stellar emission, PAH emission, 
power-law emission, and thermal grey-body emission. This is based 
on the mid-IR model used in Sajina, Lacy, and Scott (2005). 
The stellar emission is accounted for by a 10 Gyr-old, solar metal- 
licity, Salpeter IMF, single stellar population (SSP) spectral tem- 
plate generated with the PEGASE2.0 spectral synthesis code (Fioc 
et al. 1997). The specific stellar model used here is not important, 
being merely a stand-in for 'old stars'. For any population older 
than about IGyr, the near-IR spectrum looks essentially the same 
(through a couple of broadband filters), the differences being in the 
shorter wavelengths. Since our sources are largely 2 ~ 0, with a tail 
extending up to z ~ 1, this can be thought of as: all of our sources 
include some stars that are as old as the Universe at the observed 
epoch. 

We use the neutral PAH emission template of Lee & Draine (2003). 
The only modification we make is the addition of the newly discov- 
ered 17 /xm feature (Smith et al. 2004). Our broadband data do not 
allow us to investigate the probable variations in relative PAH fea- 



ture strength, and therefore any reasonable PAH template can be 
used here, with the understanding that it is only an approximia- 
tion for any given galaxy. A power-law, cx i^^", is a proxy for 
the warm, small grain emission (roughly < 60 /^m). This is cut-off 
as exp(— 0.17 x 10^* ¥Lz/u) in order not to interfere with the far- 
IR/sub-mm wavelength emission, which is described by a thermal, 
(xu''^^'^[exp{hu/kT) — component. Since we typically 

only have the 24 fim point to constrain this component, we fix a 
at 3.0, which smoothly connects this to the greybody component. 
Extinction, parameterized by rv, is applied using the 7?v = 3.1 
Milky Way-type extinction curve of Draine (2003). This includes 
the 9.7 fim Si absorption feature, which thus becomes noticeable in 
this model at high opacities. We consider only a screen geometry, 
i.e. = Jo exp(— r^). Fig. shows the above phenomenological 
break-up of the SED. We solve for the best-fit model, and the as- 
sociated uncertainties in the parameters via Markov Chain Monte 
Carlo (MCMC). Details of the fitting procedure and error analysis 
are given in Appendix A. 

This model is complicated by the lack of spectroscopic redshifts for 
about 1/3 of the sample. We discuss our approach to determining 
photometric redshifts for these sources in Appendix B. 



4 FIT RESULTS 

The quantitative conclusions of our model-fitting are presented in 
Table |2| where the errors represent marginalized 68% confidence 
limits. A quality flag is also given to indicate whether or not prob- 
lems exist with the fit. Below we discuss both the general aspects of 
the fitted SEDs, and details of the parameters presented in Table|2| 

4.1 General trends in the SEDs 

Fig. |4| shows the data overlaid with the best-fit SED model for 
sources ranging between z ~ (top row) and 2 ~ 1 (bottom row). 
As expected, there is a greater range in SED shapes than accounted 
for by the data uncertainties. However, there are natural groupings, 
for example one can easily distinguish ULIGs (e.g. Nl-064, Nl- 
078) from LIGs (Nl-015, Nl-039) and from normal galaxies (e.g. 
Nl-002, Nl-009). The grey shading in Fig.|4]gives a sense of the 
uncertainty of the SED fits. Without prior assumptions on the shape 
of the SED (s.a. by using templates) as in this case, if we do not 
have constraints on both sides of the thermal peak (e.g. Nl-078), 
a wide range of SED shapes and consequently temperatures, lumi- 
nosities etc. are acceptable. 

In Fig. |5| we address the question of what is the 'typical' 
normal/cold galaxy spectrum based on our model fits. In order to 
minimize redshift bias, we construct a composite spectrum from 
all sources with redshifts below 0.1, and for which 24 fim data 
are available. To minimize luminosity effects, we also normalize 
the spectra at 4.5 pm (which point resulted in the least amount of 
scatter across the SED) and effectively is a normalization in stel- 
lar mass. This procedure is consistent with results for ISO Key 
Project normal galaxies, where the stellar-l-PAH SEDs were found 
to be fairly constant (Lu et al. 2003). Comparing the results with a 
number of common SED models, we find that our sample is best 
described by fairly cold spectra with low mid-IR/far-IR ratios. In 
Fig.|4] we see that the LIGs (i.e. presumably starbursts) also have 
fairly cold spectra, unlike that of the 'prototypical' starburst M82 
(this was already noted in D05). 

Even more surprisingly, this cold trend does not appear to re- 
verse for the highest luminosity sources including ULIGs. It must 
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Figure 4. The SED fits for our sample, where we show both the data and best-fit model together with an uncertainty band. Note that the least well-defined 
SED is that of Nl-078, which is the only source without data in the range 24 /im- 170 fim. The sources are ordered so that redshift increases from top to 
bottom. Note the characteristic increase in the dust-to-stars emission ratio as one goes to higher redshift (and hence more luminous) sources. 



be noted however that it is not clear how common are such SEDs 
for ULIGs in general as we only have two spectroscopic ULIGs in 
our sample and our far-IR selection likely biases us toward colder 
sources. In Fig.|5|we show the best-fit SE Ds for the two spectro - 
scopically confirmed 2 ~ 0.5 -1.0 sources tehapman et alj|2002h . 
The dotted line is an Arp220 template iSilva et al. 1998) at the ap- 
propriate redshift. In both cases the sub-mm data are well fit by the 
Arp220 template, but the mid-IR data differ by an order of mag- 
nitude. Our photometric highest-z source (Nl-78) shows similar 
trends. Note that the PAH features in the best-fit Nl-064 spectrum 
are not constrained by any data point and therefore a model with no 
PAH emission is quite acceptable as well (see the spread in Fig.|4j. 
Despite it being poorly constrained for the few highest-z sources, 
we find clear evidence for prominent PAH emission for nearly all 
sources in our sample, regardless of redshift (and hence luminos- 
ity). 

The weak mid-IR continua, and strong PAH emission of our 
sources suggest that they are star-formation rather than AGN dom- 
inated (see e.g. Sajina, Lacy, Scott 2005). 



4.2 Luminosity and SFR 

In Table |2| we present both the total infrared power output of our 
sources (L^-iooo), as well as the luminosity due to PAH emission 
alone. 

The overall luminosities we derive are typically a 
few X 10^" Lq for our lower-2 targets, consistent with roughly 
L* galaxies. Lpah is obtained by integrating under the PAH 
component alone. We find the LpAn/Lm. fraction to be typically 
5 - 15%. This is consisten t with the results for the ISO Key Project 
galaxies <Dale et alliooH) . 

Since such well-sampled SEDs are rare, to obtain the overall 
infrared luminosity, extrapolations from the mid-IR are common. 
However, the coldness of our sources means that, if such sources 
are the norm, prior relations might not be applicable. In Fig.0we 
compare our results on the mid-IR/total-IR relation with previous, 
IRAS-based, results (e.g. Takeuchi et al. 2005). The solid line is the 
best-fit for our sample, which is: 

logL24 = (1.13 ± 0.05) X logLiR - (2.5 ± 0.5). (1) 
Our rest-frame 24 pm fluxes are below those expected from the 
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Figure 5. The shaded region represents the range of spectra, where we have 
included all sources with z <0. 1 and which have 24 [ixa data (except Nl- 
009, which has an unusually strong stellar component, for our sample). All 
SEDs ai'e normalized at 4.5 ^m, which is a neutral point between pure stel- 
lar and PAH emission. The shaded region should be regarded as a composite 
spectrum, representative of the 'cold' FIRBACK sources. For comparison 
we overlay a number of templates (see legend). 
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Figure 6. The best-fit SEDs for the two spectroscopically-confirmed 
higher-z sources. The IKAS upper limits are also shown, although they are 
not included in t he fit. The dotted line shows the appropriately redshifted 
Arp220 template <Silvaetalll998i) . The Spitzer fluxes suggest cooler and 
less luminous sources than previously assumed. 
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Figure 7. The 24 /.tm luminosity vs. total infrared luminosity. Note that 
within the uncertainties we find no difference between the IKAS 25 /^m 
luminosity and the MIPS 24 /im luminosity. See Section l4!2l for details. 



IKAS relation by ~ 0.4 dex. Although we have too few sources at 
the high-L end to claim this conclusively, it appears that the rela- 
tion drops even further for these sources. More likely, however, this 
is just an indication of the underlying scatter in the relationship due 
to variations in the SED shape. 

This discrepancy with the earlier relation is most likely a se- 
lection bias. The sample used for the IRAS relation is flux-limited 
to all four IKAS bands, leading to a bias toward sources with 
stronger warm continuum, unlike our 170 /^m selection where the 
bias is toward the presence of cold dust instead. We return to this 
point in Section 5.1. 

About 70% of the sample have log(I//LQ) < 11, while the re- 
maining 30% have log(L/L0) > 1 1. This can be thought of as the 
break-up between fairly quiescent and more actively starforming 
galaxies (we address the mode of star-formation of our galaxies in 
detail in Section 5.3). The total infrared lumino sity is well known 
to trace the curren t SFR of a galaxy (see e.g. lKennicut3ll998bt 
iKewlev et alj2002h . We use the Kennicutt relation here which is: 



SFR 



= (1.8 X 10"'")- 



(2) 



Moyr-i ^ ' Lq ' 

Note the slightly different definition of Lir, which is accounted 
for here (i.e. 8 - 1000 /im rather than 3 - 1000 [im). Typically our 
galaxies have SFR of ^^5 M^lyT, which is slightly enhanced with 
respect to local quiescent spirals. As expected, the SFRs rise to a 
few hundred MQlyr for our few ULIGs. 



4.3 Stellar mass and dust obscuration 

As part of our SED model fitting we estimate the fc-corrected, and 
dust corrected near-IR stellar emission. The near-IR is preferred for 
deriving stellar mass as it is fairly robust against both uncertain- 
ties in the dust obscuration level (a factor of ~ 10 less so than in 
the optical), and to the details of the stellar population and SFH 
(being largely sensitive to the red giant population only). With- 
out sufficient optical coverage, we cannot properly account for the 
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possible contribution of young stars. We merely calculate the un- 
extincted, restframe Lk, and convert this to stellar mass by assum- 
ing an appropriate A'-band mass-to-light ratio, ')k- Despite its rel- 
ative insensitivity to young stars, the near-lR mass-to-light ratio 
does increase as the fraction of young/massive stars in the galaxy 
increases. In addition, a degeneracy exists in that a more massive 
and dustier galaxy can appear similar to a less massive and less 
obscured one. The significance of this effect grows for sources 
without y-band data or with poor near-IR SNR, both of which are 
more common for the higher-z sources. However, since the optical 
depth is a free parameter in our model, it is already accounted for in 
the resulting Markov chains. The uncertainty in mass-to-light ratio 
however is not included. The additional uncertainty on the stellar 
m ass is: d log Af, = d log 7. A mean value of '^k = 0.6 was found 
b^jBel^^^^onaj200l|) for their sample of spiral galaxies, while 
lGi^^a^^Ln2000r fin d 7x to be ^^0.92 fo r their sample of 
starburst and Hll galaxies. iBell & de Jonj ('2001") consider a wide 
range of reasonable population models, finding d log jk ~ 0.2 dex 
(although this is smaller if the B-K colour is known, due to the 
effect of the young stars). The SFRs derived for our sources in 
the previous section suggest enhanced star-formation activity com- 
pared with normal spirals. Therefore, the Gil de Paz et al value 
of jK = 0.92 is likely more appropriate for our sample, while we 
assume the Bell & de Jong uncertainty (from the unknown SFHs) 
of 0.2. Adding this in quadrature to the MCMC derived errors 
suggests that more realistic uncertainties here are about twice the 
quoted ones. 

With the above set of assumptions, we find 
{A//star)~2x 1O"M0,2 with the log(L/L0)<ll sources 
and those with log(I//L0)>ll differing by ~0.6dex. Note, 
however, that from the above discussion we expect that somewhat 
higher values for 7 are appropriate for the more luminous (higher 
SFR) galaxies, and vice versa for the lower luminosity ones. 
Taking this into account, the observed difference in mass may be 
even largely in reality. 

For most sources, we find modest levels of dust extinction 
(0 < TV < 5, peaking at ^ 3) for most sources. This assumes a 
screen geometry, and would increase in a uniform mixture of dust 
and stars. The amount of extinction generally increases for 2: ^ 0.3 
sources, with N 1-040 requiring the greatest optical depth, con- 
sistent with its near-IR faintness (note however that the mean- 
likelihood and MCMC approaches disagree on the optical depths 
of the high- 2: sources - see Appendix A). 

The values of tv we find for the bulk of the FIRBACK Nl 
galaxies, is consistent with the average rv ~ 3 found from the 
Ha/Hf) ratio of the FIRBACK-South galaxies (Patris et al. 2003). 
This means that, typically, the optical depths to which the light of 
old and young stars is subjected do not differ dramatically when 
integrated across the galaxy (see also Section 5.2). 



4.4 Dust properties: T - /3 relation and dust mass 

The single greybody approach we take (see Section IT2t in mod- 
eling the dust emission of these galaxies, although in common use 
and the only one possible when just a few data points are avail- 
able, is clearly a n approximation (see also lDunne & Eale3l200lt 
iBlain et aljboOSh . More realistically, a distribution of dust grain 
characteristics, such as optical properties and geometry, results in 



^ For comparison, typical M* values are ~ 3 X 10^'' Mq. 
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Figure 8. The best fit va lues (black dots) for Td and /3. The solid line 
is the iPupac et al J l2003l) relation (see Appendix C). For comparison the 
grey scale points show the X^in ■'^gion for N 1-024. The error bars rep- 
resent the average 68% uncertainties on individual sources (note that it is 
driven by the least constrained sources, while the grey points give a better 
sense of the scatter in well constrained sources). The outlier near 70 K is 
Nl-078, for which no data exist between 24 fim and 170 /xm. 



a distribution of effective temperatures and /3's. In addition, degen- 
eracies arise in fitting this model due to the functional form and 
spectral sampling (see Appendix Iclf or further discussion). 

For our sample, values of Td ~ 20 - 30 K and high values of 
P (~ 2) are most representative. These are in good agreement with 
the typical far-IR temperatures (~ 22 K) of local spiral and irreg- 
ular galaxies 1 C ontursi,.2001.; Stickel et alii200 0). as well as the- 
oretical expectations which place the big grain emission (domi- 
nating the thermal peak of normal galaxies) of a standard ISM at 
Td ~ 15 - 30 K with /3 ~ 2 ( Draine & Lee 1984). Our estimates are 
also consistent w ith, but somewhat cooler than, the values found by 
iTavlor et alJ i2005»^ for the FIRBACK galaxies in ELAIS-N2. 

The cold dust mass is esti mated in the usual way: Md = 5*850 * 
Dl[{l + z)HiB{T, i/e)]"^ (e.g. lFarrah et alj2 002'). where B(T, Ve) 
is the Planck function for the given temperature, at the emis- 
sion frequency, and the dust absorption coefficient k parame- 
terizes the unknown grain properties. We as sume the valu e of 
K=0.077±0.02m^kg"\ as advocated by Jam es et alj j2003h . al- 
though there are arguments for a value up to 3 times higher 
iDasyra et al. 2005J. We return to the effect of increasing k in Sec- 
tion 5.2. We find Md ~ lO'^-lO* M©, which is comparable with 
the masses found for the SCUBA Local Univ erse Galaxies Sur- 
vey (SLUGS) galaxies (Dunne &Eales 2 00 ll) . As expected, the 
higher z, more luminous galaxies have higher dust masses, typi- 
cally - 10** - 10^ Mq. 



4.5 Size and morphology 

The closest galaxies in our sample appear disk-like by visual in- 
spection of the IRAC images (see Fig. 0. In addition, we would 
like to determine the typical sizes of our sample, and whether there 
is a noticeable morphological difference as a function of infrared 
activity (i.e. Lir). Fortunately, all 22 sources considered here are 
at least mildly resolved by IRAC (diffraction limit ~ 2" at 8 ^m). 
The effective radii given in Table |2| are the 8 /im half-light radii 
converted to physical distance (note there is no noticeable differ- 
ence here if we used any other IRAC band). We find that, typically, 
^?eff ~ 5- lOkpc, although a tail of the sample extends to larger 
radii. These sizes are consistent with those of large spirals. 
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5 ANALYSIS 

In this section, we use the results of the previous section in at- 
tempting to arrive at a consistent picture of the physical proper- 
ties of our sources. We particularly focus on three related ques- 
tions. Why are all our galaxies, including the ULIGs, colder than 
expected? Can we say something about the spatial distribution of 
the star-formation activity - i.e. is there indication of centrally con- 
centrated starbursts surrounded by older stellar haloes, or is the 
star-formation by contrast distributed throughout the galaxy? And 
finally, is the typical FIRBACK galaxy actually starbursting?. 



5.1 Why cold LIGs and ULIGs? 

Both our general inspection of the SEDs and the derived dust tem- 
peratures for our sample suggest that fairly cold spectra with weak 
(compared with IRAS galaxies) mid-IR emission describe our en- 
tire sample, despite three orders of magnitude variation in bolomet- 
ric luminosity. This can be seen directly by examining the L - Td 
relation (see B lain et al, 2003) . This relation is of particular sig- 
nificance to far-IR/sub-mm selected samples, as in principle it can 
reveal something of the nature of the sources based on only a few 
spectral points (enough, for example, for a single greybody fit), be- 
cause it is simply a relation between the location of the far-IR peak 
and its overall strength. In Fig.|9| we compare the location of our 
sample in the (L, Td) plane with other infr ared-selected samples , 
namely the SCUBA-selected high-z sources tehaoman et al."2005^ 
and the bright 7RAS-selected galaxies (Dunne & Bales 2001). We 
overplot empirical relati ons appropriate for merging and quiescent 
galaxies <Bamardl200l . These can be understood intuitively by 
recalling that ideally, L cx R'^T'^. This leads to luminosity increas- 
ing with temperature if the size is kept constant, or conversely to a 
more concentrated starburst being hotter than a galaxy of the same 
luminosity but with more spread-out star-formation (as in through- 
out the full disk). With a few exceptions, our sample is consistent 
(within the uncertainties, and given the different selection from that 
of IRAS galaxies) with the relation for quiescent galaxies (for fur- 
ther discussion of the relation see C hapman et al. 2003) . and con- 
sistent with our size results in Section 4.5. 

The above discussion, taken a step further, suggests 
that our cold ULIGs might show extended star-formation 
activity, rather than the usual nucle ar-concentration see n 
in the late stages of major mergers iVeilleux et alJ l2002h . 
Such a scenario has in particular been proposed f or some 
SCUBA galaxies I'Kaviani. Ha ehnelt. & KauffmaniJ l2003t 
[sfstathiou & Rowan-Robinson 2003). Locally, only nuclear star- 
bursts appear capable of reaching the requisite SFR/luminosity. 
However, at high-z, where presumably the progenitors are less 
evolved and more gas rich (i.e. with lower stability thresholds), it 
is possible that a major merger results in a disk-wide starburst (or 
rather many pockets throughout the disk) instead (Mihos 1999). 
Alternatively, a cause of cold ULIGs might be extreme opacity 
around a nuclear starburst and/or AGN. 



5.2 The spatial scales of the old and young stars 

In the previous section, we argue that the coldness of our sources, 
including both LIGs and ULIGs, might be due to star-formation 
spread over an extended cold disk, rather than a nuclear starburst 
surrounded by an old stellar halo. Here we apply two independent 
tests of this scenario. 

Our first test is based on the assumption that the IR emission is 
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Figure 9. A comparison of the luminosity-temperature relation for our 
sample (solid cir cles) with other infrared -selected galaxies including IRAS- 
bright g alaxies iDunne & EaleJ l200ll) {crosses), and SCUBA-selected 
galaxies iChaDmai^^^200^r fopew squares). The solid and dashed lines 
repres ent the loci for merging and quiescent galaxies respectively i Bama^ 
I2OO2I) . Note that our sample is consistent with the expectation for quiescent 
starformers. 



dominated by power absorbed by the dust from young stars, while 
the dust responsible for the near-IR opacity traces the old stellar 
population. The two should match in a disk-wide star-formation 
picture, while there is no reason for a correlation in a nuclear star- 
burst. In Section 4.4 we derived the dust mass based on the far-IR 
emission of our galaxies. In principle, the optical depth derived in 
Section 4.3 can also be converted to dust mass, provided the sur- 
face area of the optical region is known. We do this usin g the fol- 
lowing equation, adapted from lLisenfeld & Ferraralil998h by using 
rs = 1.3 TV and assuming the geometry is an inclined disk: 

M°P* = (1.1 X 10^)TvR^cosi, (3) 

where i is the the inclination, such that 0° is face-on and 90° is 
edge-on. 

In Section 4.5, we presented the effective radii of our sources 
as derived from the IRAC images. We found no difference between 
the 4 IRAC bands and therefore it is safe to assume that these radii 
represent the spatial extent of the old stellar population. In Fig. 110! 
we compare this optically derived dust mass and the IR-derived 
dust masses. The two clearly correlate, as expected. However, the 
best-fit slope is 0.86 unlike the expected 1.00. The arrows define 
the directions in which either increasing disk inclination or increas- 
ing K push the points. Note that increasing inclination is equivalent 
to decreasing the size of the region contributing the bulk of the 
IR emission (i.e. towards nuclear starburst). Therefore at the high- 
mass end, the good agreement suggests that the dust (behind the 
IR emission in particular) is distributed throughout the entire disk 
rather than being concentrated in the nuclear region (see next sec- 
tion). Somewhat contradictory, the high values of tv we find are 
also fully consistent (but the lower values implied by the likelihood 
analysis are not - see Appendix A). The two can be reconciled by 
a more complex geometry. 

At the low mass end we find consistently higher dust masses 
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Figure 10. Optically-derived dust mass (from absorption) compared with 
IR-derived values (from emission). The solid line is the best-fit, while the 
dashed line is the expected y = x line (for face-on disks). The directions of 
increasing disk inclination, and k are shown. 



than expected from the above formula and the derived optical 
depths (hence the shallower slope). This may merely indicate that 
the highest-SFR regions responsible for the dust emission are be- 
hind a foreground screen of much less obscured older stars which 
are responsible for the near-IR emission. Another interpretation, is 
that our IR-based dust mass is ov erestimated. Ad opting the extreme 
end of K values in the literature iPasvra et alJ20 0.5\ k~0.2, will 
decrease our estimates by about 0.4 dex, making the two estimates 
much more consistent. The inclination cannot reconcile them as it 
only affects the points to the right of the dashed line (face-on case). 
This relation should probably not be over interpreted, since both ex- 
pressions are approximations to much more complex systems and 
seeking a perfect agreement between them is therefore unrealistic. 
Within the uncertainties, we are gratified that our two independent 
estimates of the dust content agree as well as they do. 

Given the above uncertainties, and especially if k is underesti- 
mated, it is still possible that some of our sources (especially among 
the higher-L ones) are primarily powered by nuclear starbursts. We 
can determine the actual sizes of the star-forming/bursting regions 
in our galaxies by using the well kno wn relation bet ween gas sur- 
face density and SFR surface density JSchmidJl959l) . The implied 
universality of the star-formation efficiency has been demonstrated 
observationally over many orders of magnitude in SFR density 
iKennicutt.l998bll) . The slope of the relati on SFR ocp" is n= 1.4 
(the Kennicutt relation: |Kennicutlll998d) . However, the depen- 
dence of both quantities on radius is the same (since SFR and gas 
mass are derived from integrated properties, their respective sur- 
face densities are simply oc 1/r^), meaning that a change in the 
effective radius takes the form of a translation along a line of slope 
n = 1.0. Thus, assuming that the efficiency of star-formation is in- 
deed universal, any departure from the Kennicutt relation could be 
attributed to a difference in effective radius. 

We estimate the gas mass by assum ing the Milky Way gas-to- 
dust ratio of 100 iKnapp & KerJll974h . although the derived val- 
ues, for external galaxies especially, show a large spread of about 




Figure 11. The Schmidt law, Sgpa = a Sgj^j,, where observationally 
n= 1.4 lKennicut3ll998hl solid line). Our galaxies are plotted as the filled 
circles. The errorbars represent their average 1 a uncertainties. For compar- 
ison we overplot the Kennicutt spiral galaxies (open triangles), and starburst 
galaxies (open stars) samples. The inset shows the size distribution for our 
sample. The slope = 1 arrow shows the direction of increasing radius, while 
the horizontal arrow shows the direction of increasing dust-to-gas ratio. 

an order of magnitude (see e.g.'Stic kel et alllOOd) . With these es- 
timates, in Fig. 111! we overplot our sources onto the original Ken- 
nicutt sample of normal spirals, and IR-bright starbursts. Most of 
our galaxies lie on this relation (within the errors) somewhat above 
the spirals, consistent with the SFRs found in Section 4.2. This sug- 
gests that the bulk of our sample is indeed forming stars throughout 
their extended disk, although slightly more actively (due to more 
gas-rich disks) than in local spirals. 

The good agreement between our sample and the Kennicutt rela- 
tion suggests that both the estimated radii and the assumed standard 
dust-to-gas ratio are approximately correct. The inset in Fig. 1111 
gives the size distribution of the sample, which is roughly consis- 
tent with that of local spiral galaxies ( Vertchenko & Ouiroga 1994)- 
An HST-based study of the size distributions of slightly higher red- 
shift (z ~ 0.2- 1) disk galaxies suggests effective radii in the same 
range with the peak moving from ^ 6 kpc to ~ 4 kpc as the redshift 
decreases iRavindranath et all2004b . 

As a final test, we double check the above conclusions by 
looking at the relation between Lm, and T'^'^^ (which should 
scale with luminosity, for simple geometries). As expected we find 
a good correlation between the two. The linear fit gives: 

logL = (0.9 ± 0.1) X log(T*+''j?^) + (1.5 ± 1.3), (4) 

where L is in Lq, T is in K, and R is in kpc. There are no obvi- 
ous outliers (the rms is 0.4), again suggesting that given the cool 
temperatures derived, the radii inferred from the 8 images are 
appropriate (to within a factor of 2). 

5.3 Starbursts? It's about timescales 

In the previous section, we found that the typical FIRBACK galaxy 
has modest SFR, typically ~ 5 M0yr~^, but this is somewhat en- 
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Figure 12. Stellar mass vs. specific SFR. This compilation of compari- 
son galaxies is from lGildePazetalJ i200(D : howev er we have conver ted 
their SFRs to the standard Kennicutt et al. relation iKennicutll [19983) (a 
difference of -0.87 dex). Shown are: normal spirals (open circles); local 
Ho-selected star-forming galaxies (crosses); dwarf/HlI galaxies (stars). 
Our FIRBACK sample is shown as the filled circles. Note that for con- 
sistency, we have co nverted our IR-based SFR to Ho-based SFR using the 
iKewlev et alJ 12002 ) relation. The aiTow shows the direction of increasing 
Ha luminosity from Gil de Paz et al. (2000). The mass boundary (shown 
by the vertical line) is the observed bou ndary in physical properties seen in 
the SDSS data iKauffmann et all2003l) . 



hanced activity compared with local quiescent spiral galaxies. Our 
investigations so far suggest disk-like star-formation, with no indi- 
cation of mergers, especially for the low- 2; sample. Now we would 
like to address the question of whether this enhanced star-formation 
could reasonably be described as a 'starburst'. There is no single 
absolu te definition of a starburst, but some commonly used indi- 
cators taeckman 2005^ include: a high intensity of star-formation; 
a high birthrate parameter (SFR/(SFR)); or a low ratio of the gas 
depletion timescale to the dynamical timescale. The question is - 
when is 'high' high enough, and when is 'low' low enough. We 
now discuss each criterion in turn. 

The intensity of star-formation, was indirectly addressed in 
Fie. 111! Our sources appear to fall somewhere in-between quies- 
cent spirals and nuclear starbursts (the Kennicutt et al. starburst 
sample are all nuclear). Here the highest-z sources (Nl-064, and 
Nl-078) are apparently the only unambiguous starbursts. The re- 
gion occupied by the bulk of our entire galax ies is in fact a lso the 
locus of the centres of disk galaxies alone IlKennicuti 1 998bf) . So, is 
the star-formation in our galaxies fairly smoothly distributed across 
the disk or in a collection of bursting knots sprinkled throughout? 
Are we diluting our results by assuming a smooth disk here? It ap- 
pears that the intensity argument is inconclusive when applied to 
unresolved galaxies except in the most extreme cases. 

The birthrate parameter is slightly misleading, because of its 
dependence on the epoch of observation (although this is not too 
severe for our sample). A r elated quantity, which w e examine next, 
is the specific SFR (see e.g. lGil de Paz et all200Cl) . which is equiv- 
alent to the birthrate assuming a uniform epoch. The best way to in- 
terpret this is through comparison with other samples whose prop- 
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Figure 13. The ratio of gas depletion to dynamical timescales. The two 
histograms are for the case where dark matter is not included (solid line) 
or when it is (dashed line). The shaded region corresponds to the classical 
definition of a starburst. 



erties are already known, and for this purpose we use the compila- 
tion of lGil de Paz et alJ <2000h . However, such comparisons should 
be carried out with some caution, since the samples are selected 
in very different ways and their quoted SFRs are based on differ- 
ent indicators. The main distinction here is that these other samples 
are based on the strength of the Ha emission line, while ours are 
based on the overall IR emission. The two are comparable, but not 
identical (largely due to the effects of dust extinction on the line 
strength at the high-L end and possible contribution of older stars 
to the IR e mission at low lumin osities). Here we use the empirical 
relation of iKewlev et alj i2002h to convert our IR-based SFRs to 
Ha-based ones. Lastly , we standardize the Gil de Paz et al. SFRs to 
the Kennicutt relation ( Kennicutt 1998a). In Fig. lI2l we plot stellar 
mass vs. specific SFR for our sample. For comparison we over- 
plot the data for lo cal normal spirals, stron g Ha emitters, and Hll 
dwarfs, taken from lGil de Paz et alJ i200ol) . Fig. 1 121 suggests that, 
overall, our sample is somewhat less massive and more active than 
typical local spiral galaxies. However, the bulk are larger and less 
star-bursting than the average local Hll dwarf. In fact they most 
strongly resemble the comparison samp le of Ha-selected local star- 
forming galaxies iGil de Paz et alj200 0). This is to be expected, as 
the available spectra of our sources all show prominent emission 
lines (D05). The two outliers are the two ULIGs Nl-064, and Nl- 
078 which are predictably both massive and active. 

We finally turn to the ratio of the gas depletion timescale 
(Tdcpi) to the dynamical timescale (rdyn) of a galaxy. Unlike the 
above two considerations, the starburst criterion is clearer here. If 
the current rate of SFR is such that, if sustained, the available gas 
reservoir will be exhausted in a time shorter than the dynamical 
time, then the source is bursting. The depletion timescale is sim- 
ply Mgas/SFR. We find that, typically, Tdcpi ~ 8x 10® yr. The dy- 
namical time is defined as rdyn ^ {R/GA4tot)^^^ ■ To estimate 
this, we begin by neglecting the dark matter contribution, and cal- 
culate the dynamical time using the gas mass plus the stellar mass. 
We then note that including any dark matter contribution would 
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only have the effect of decreasing Tdyn and thus any source that 
appears quiescent will remain so, regardless of the extra mass (we 
return to this point below). Using this approach we find that typ- 
ically Tdyn ~ 1 X 10* yr. Fig. E| shows a histogram of rdcpi/rdyn 
using estimates both with and without dark matter in Aftot . When 
the dark matter halo mass is included^, the only starbursts in our 
sample are Nl-013, and Nl-078 (with Nl-064, and Nl-077 being 
borderline). Throughout, we assume our IRAC-based sizes, which 
we argued are the effective scales of the star-formation activity (see 
Section 5.2). The bulk of our sample, despite their being apparently 
somewhat more active than local, quiescent spirals, do not appear 
to be bona fide starbursts. 



6 DISCUSSION 

In this paper, we have presented a comprehensive study of the full 
infrared SEDs of a sample of 22 FIRBACK galaxies. Our new 
multi-component model together with the MCMC fitting technique 
have allowed us to test various previous assumptions about FIR- 
BACK galaxies, as well as the spectral templates used to model 
them and related populations. 

We found that our galaxies have cool SEDs (low mid-IR/far- 
IR ratios), which remains true even for the higher luminosity 
sources (LIGs and ULIGs). It is not yet clear how common such 
sources are, and how much of a concern are they for galaxy 

evolution models. However, the Spitzer 24 am number counts 
I I I 

iPapovich et al . 2004) already suggest that the existing SED spec- 
tra might not be accurate in the mid-IR for the z ~ 0.5-2.5 lumi- 
nous and dusty sources thought to comprise the bulk of the CIB. To 
correct for this, the necessity for 'downsizing' the mid-IR in their 
starburst spectra was already pointed out by Laeache et al. ( 2004). 

One of the primary questions we wanted to answer with this 
study was the nature of the brightest galaxies contributing to the 
CIB. Initially from the sub-mm data (S03) and later through optical 
spectroscopy follow-up (D05), it was already known that the bulk 
of the FIRBACK galaxies are \ow-z, moderate luminosity sources, 
with a small fraction of 2 ~ 1 ULIGs . Beyond that however, lit- 
tle has been known of their nature, and in particular their masses, 
sizes, and modes of star-formation were poorly constrained. Here 
we have addressed these issues from a number of (admittedly not 
always independent) perspectives and concluded that the bulk of 
the sample is consistent with ~ M* mass galaxies, which are form- 
ing stars somewhat more actively than local spirals, but are likely 
not actually starbursting according to the usual definition. Their 
mode of star-formation is consistent with slightly enhanced activ- 
ity in the disk, rather than the nuclear bursts associated with ma- 
jor mergers. Their cool colours also exclude any significant con- 
tribution fr om AGN acti vity. A study of the ELAIS-N2 FIRBACK 
galaxies bv lTavlor et all ^005) used a set of theoretical templates 
for cirrus (i.e. quiescent), starburst, and AGN galaxies. They con- 
clude that 80% of the FIRBACK sources are starbursts, while 20% 
are cirrus galaxies. Since we find enhanced star-formation, but use 
more stringent starburst criteria, our conclusions are consistent. 



3 Using the relatio n A/haio/Mgas 35(Afgas/lO'^M0)-' 
iMacLow & FerrarJl998ft 




1 2 3 4 5 
redshift 

Figure 14. Detectability of FIRBACK-like galaxies as a function of 
redshift, compared with the Spi tzer 24 /jm surve y limits for SWIRE 
ISurace et all2004 . and GOODS iCharv et all2004 . The dashed line cor- 
responds to the SED of a typical low-z FIRBACK source (Nl-007), while 
the solid line represents a LIG (Nl-015), and the dot-dashed line is a cold 
ULIG (Nl-064). Note that the level of PAU contribution is uncertain in the 
last case and therefore it might not be detectable by SWIRE at z ~ 2 as 
suggested in this figure. 

6.1 Implications for the faint 24 /im sources 

Over the past year, Spitzer has revealed large numbers of faint 
(roughly ~ 10-1000 /ily) 24 /^m sources. There has been much 
speculation about their nature, with obvious implications for var- 
ious galaxy evolution models. In Fie. 1141 we evolve the SEDs of 
an L* galaxy, a LIG and an ULIG from our FIRBACK sample. 
We compare these with the detectability thresholds of the wide and 
shallow SWIRE survey 1 Surace et al. 2004) a nd the deepest current 
Spitzer survey, GOODS-North iCharv et aL"200?). The SWIRE 
sources are expected to be predominantly low-2:, meaning that 
FIRBACK-like galaxies (at only slightly higher redshifts) make up 
a significant fraction of them. For the deeper GOODS survey, FIR- 
BACK analogues are expected to be a large component up to 2 ~ 1, 
while in the 2^1-2 range GOODS appears to be dominated by 
11 < log{L/ Lq) < 12 sources, such as the ULIG-tail of our sam- 
ple, as discussed bv lCharv et alJ i2004h . Overall, Fig. ^] suggests 
that understanding the faint 24 /im sources requires taking into ac- 
count not only M82 or Arp220 type sources, but also the much less 
luminous and colder FIRBACK sources, especially in the crucial 
2 ^ 1.5 regime. This is supported by the fact that while SCUBA- 
selected sources have close to 100% detectability in the cur- 
rent generation of 24 /im surveys fe gami et alj|2004i iFraver et alJ 
2004; Pope et al. 2006), most of the 24 /xm- selected sources do 
not have individual SCUBA-counterparts (stacking a nalyses show 
them to have individual 850 /im flux ~0.5mly 1 Serjeant et alJ 
2004,) . which is below the confusion limit of SCUBA. 



6.2 Implications for the SCUBA galaxies 

Fig. 1141 indicates that if the SCUBA galaxies are as cold as our 
ULIGs (and given that their redshifts are typically in the range 
2 ~ 2 - 3) then relatively shallow Spitzer surveys (e.g. SWIRE) 
cannot detect any but the most luminous or the most AGN- 
dominated ones among them. This has implications for the Spitzer 
overlap with SCUBA galaxies in w ide-field sub-mm surveys, such 
as SHADES ^ Mortier et ai]l2005ll . and those being planned with 
SCUBA2. The difference between the FIRBACK far-IR and blank- 
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sky sub-mm sources is further highlighted in Fig llSI Here we show 
shaded regions to indicate the colours of normal log L < 1 1 galax- 
ies, as distinct from log L > 1 1 ones, indicating their rough red- 
shift evolution. The influence of the PAH feature complex as it 
traverses the 24 /im filter at z ~ 2 is evident. It appears that, de- 
spite some overlap with our cold ULIG template, the bulk of the 
SCUBA sources are redder in both colours. As the 24 /im flux pulls 
them in different directions here, this rather suggests that the red- 
dening of one colour drives the overall effect. The S24 / ^8 colour 
can be reddened by redshift, optical depth, the presence of AGN, 
or increased PAH strength (at 2 ~ 2). Fully addressing the SEDs of 
SCUBA galaxies is clearly beyond the scope of this paper (see Pope 
et al. 2006). Here we merely show one scenario, i.e. that increased 
level of obscuration can account for the SCUBA galaxies' colours. 
We take our best-fit Nl-064 SED (tv = 5), and subject it to addi- 
tional extinction to a total of tv =14. This is shown as the dashed 
curve in Fig. 1151 The most SCUBA galaxies now fall in-between 
the original and revised Nl-064 template; however, we remind that 
we have not fully explored the parameter space, and this is merely 
a consistency argument. 

The three FIRBACK ULIGs we have studied here are Nl-040, Nl- 
064 and Nl-078. In Section 5.1, we concluded that these sources 
are colder, and likely have more extended star-formation activ- 
ity than typically assumed for ULIGs. This is qualitatively con- 
sistent with the cirrus models of Efs tathiou & Rowan-Robinson 
iEfs tathiou & Rowan-Robinsoni [2003h . which was also claimed as 
a possible model for the SCUBA galaxies. A small fraction of 
the SCUBA sources are indeed consistent with the SEDs of these 
galaxies, redshifted slightly to z^l-2. Most SCUBA galaxies, 
as discussed above, are consistent with the conventional view that 
they are highly obscured, extreme starformers (possibly including 
AGN) - the result of major mergers. 

It is also worth noting that the use of the proposed 850 to 
24 fim flux ratio as a redshift indicator will at best be highly unre- 
liable, due to the 24 fim band traversing the PAH emission and Si 
absorption features in the mid-IR, as well as because of the range 
of possible far-IR SEDs. On the whole the SCUBA galaxies appear 
to exhibit a large enough range in SED shapes that applying a sin- 
gle model to their mid-IR photometry would be misleading. More 
comprehensive studies of the multi-wavelength properties of larger 
samples of sub-mm selected galaxies should help us understand the 
difference between the sources comprising the sub-mm background 
and those responsible for the CIB. 



7 CONCLUSIONS 

In this paper we have combined archival Spitzer observations of 
the ELAIS-Nl field with prior near-IR, far-IR, and sub-mm data in 
order to study the full ~ 1 - 1000 fim SEDs of the brightest con- 
tributors to the CIB at its peak. A novel MCMC fitting technique 
and a phenomenological SED model are used to make optimal use 
of the available data. Below we highlight the principal results of 
our study: 

• Contrary to expectations, we have demonstrated the exis- 
tence of vigorously star-forming sources, which nevertheless have 
low mid-IR/far-IR ratios. Our sample extends over two orders of 
magnitude in luminosity and includes surprisingly cold ULIGs. 
This is likely the result of our far-IR selection. 

• We discuss some of the issues inherent in interpreting 




Figure 15. The evolution of two diagncstic colours: 5850/524 and 
S24/SS. The solid circles are our FIRBACK sample, while the open circles 
are S CUBA-selected galaxies from the GOODS-North field IPope et alj 
l200fil) . while the cro sses are other SCUBA-selected samples from shallower 
surveys jEg^^Taf 2004; Prayer et al. 2004). The dark grey region repre- 
sents log L ^ 1 1, while the light grey region shows log L ^ 1 1 (the bound- 
ai'ies should be regarded as fuzzy). The redshift evolution of the SEDs is 
carried out in steps of 0.5 from ^ = to z = 4. The SED curves used, from 
top to bottom, are: Nl-064, Nl-029 and Nl-009, and we have used them 
to trace the approximate boundaries for redshift evoltion, split into 4 bands 
as described in the legend. The dashed curve is Nl-064 with an additional 
extinction (see Section 6.2). 

empirical SED models. In particular, conclusions may well be 
influenced by the spectral sampling. Here we emphasized how 
limited spectral sampling and intrinsic model degeneracies might 
be responsible for the claimed physical Td ~ /3 correlation. 

• A number of basic parameters for our sample are derived. 
The typical stellar mass of our galaxies is fewx 10^" Mq to 
fewx 10" M0 for the handful of ULIGs. Dust massses were 
found to be IO^-IO^Mq, where we find general consistency 
between emission and absorption derived values. Typical SFRs are 
~ 5 M0/yr. We look at combinations of the above in order to 1) 
check for consistency given the many inherent assumptions, and 
2) allow for more meaningful comparisons with local optically- 
selected galaxy samples. 

• In the specific SFR vs. stellar mass relation, our sample 
appears to have enhanced SFR with respect to local spirals, but 
below dwarf Hll galaxies. The closest match to our galaxies 
appear to be Ha-selected galaxies. This is consistent with the 
emission-line spectra of our galaxies, where available (D05). 

• All our galaxies obey the Schmidt/Kennicutt relation. They 
appear to have enhanced star-formation activity with respect to lo- 
cal spirals, but less than nuclear starbursts. The more IR-luminous 
galaxies in our sample have higher star-formation intensity, as 
expected. 

• We looked at the ratio of the gas depletion timescales to 
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dynamical timescales for our sample, and find that the FIRBACK 
galaxies fall short of the traditional starburst definition where the 
gas is exhausted on timescales shorter than the dynamical time. 

• FIRBACK-like galaxies are potentially a significant compo- 
nent of shallow Spitzer surveys, whose redshift distributions peak 
at 2; < 0.3 (e.g. FLS). Our sample does not allow us to address the 
Universal role of cold ULIGs, although it appears that sub-mm 
blank-sky surveys might also be sensitive to similarly cold (ex- 
tended and/or highly obscured) sources. 
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APPENDIX A: MCMC FITTING AND ERROR 
ESTIMATES 

Al Description of our MCMC SED-fltdng method 

We us e a Markov Chain Monte Carlo (MCMC) (see e.g. lGamermanI 
Il997h approach in order to sample the posterior parameter distribu- 
tions, allowing us to find the best solutions as well as the asso- 
ciated errors on the derived parameters. Taking a to be the array 
of parameter values defining a given model, and y to be the ar- 
ray of data values available, according to Bayes theorem, we have 
p(a|y) cxp(y|a)p(a) where p(a|y) is the posterior probability dis- 
tribution, p(y|a) is the likelihood of the data for the given model, 
and p(a) is the prior. Normalizing this (via the global probability of 
the data) is usually impossible in practice, since it requires knowl- 
edge of all possible models. Here, as is often the case, a single 
model is assumed and what we are interested in is the shape of the 
distribution p(a|y), such that the most probable values for the pa- 
rameters and their errors can be derived. What we really want then, 
is a measure of the relative probability of a given parameter value 
compared with some other value. Rearranging Bayes theorem, we 
have: 



p(a|y) ^ p(y|a) p(a) ^^^^ 
P(a'|y) p(y|a') p(a')' 

The ratio of priors is equal to 1 if a flat prior is assumed or is equal 
to e"'^"» if a logarithmic prior is chosen instead. Note that eq- 
n lAll can be used to adjust the results of a chain sampled with flat 
priors to test t he effect of other choi ces of priors ('importance re- 
sampling', e.g. lLewis & Bridle 200?). This is the approach we take 
when returning to the effect of priors in Section A.2. For now, flat 
priors are assumed for all parameters, and therefore from now on 
the posterior and likelihood distributions are used synonymously. 

In principle, arbitrary shapes of the posterior distribution can 
be sampled using a simple Monte Carlo approach. However for 
multi-dimensional problems, where the ratio of high-probability 
volume to total volume is very small, this can quickly become 
computationally prohibitive. The basic idea behind MCMC is to 
effectively sample this distribution by building up chains of ran- 
dom guesses of parameter values, where each successive guess is 
chosen from some much smaller proposal distribution, q, around 
the previous chain link. This move is accepted or rejected accord- 
ing to some criterion, which both pushes the chain toward higher 
probability regions, and allows for some random deviation from 
the straight gradient descent-type path. We follow the Metr opolis- 
Hastings algorithm I Metropolis et alll953tlHastingslll97(t) where 
all parameters are varied at once, and a guess is accepted according 
to the criterion: Ofi+i = mm[u, pi+iqi+i / piQi], where the stochas- 
tic element is provided by the random number we [0, 1]. For the 
proposal distribution, we use a multivariate Gaussian, which, be- 
ing symmetric, leads to qi+i/qi of unity. Note that the exact shape 
of the proposal distribution is not important, but its effective width 
strongly affects the efficiency of the MCMC (see below). We as- 
sume that the likelihood of a given solution is the usual expression 
given by: 

log £ = const - ' , (A2) 

where the second term is the ■ The probability (p) that the sys- 
tem finds itself in a given state is given by the Boltzman factor, 
where the 'energy' is x^, and thus pi+i/pi is exp(— Ax^/T) (this 
is called the 'odds' of the given solution). The temperature, T, has 
effectively the same function as the width of the proposal distribu- 
tion, in that it determines how easy it is for the system to jump a 
particular distance from its current state. 

Since brute-force MCMC is a fairly slow procedure, rather 
than initialize the chain at some random point we begin with some 
reasonable guess at the best model. Other possibilities would in- 
clude simulated annealing, such as used in lSaiina et alj i200^ . or 
equivalently using variable widths of the proposal distribution, or 
using any other optimization technique to find the high probability 
regions quickly. For our purposes here, starting with a reasonable 
guess is deemed sufficient since we still explore the full region of 
physically plausible solutions. 

The numerical parameters which need to be set are: the width 
of the proposal distribution for each parameter (the g-width); the 
temperature; and the overall length of the chain. Too low a g-width, 
will tend to acceptance of too many trials, while too high a value 
will give some jumps far outside the high probability regions, re- 
sulting in low acceptance rates. Trial-and-error has shown that an 
acceptance ratio in the range 10-30% is reasonable (this is sup- 
ported by empirical studies which show that ~ 25% acceptance rate 
in problems with > 2 dimensions minim izes internal correlations 
in the resulting chain jRobertser^Il993)). To find the appropriate 
g-width, we make shorter (30,000) runs, where we vary only one 
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parameter at a time. We start with a guess at the g-widths for each 
and incrementally adjust them until the resulting acceptance rates 
are 50-80% (which is appropriate for 1 -dimensional problems). By 
trial and error we have found that this leads to acceptance rates in 
the desired 10-30% range when all seven parameters are varied si- 
multaneously. This approach works better for our sample than fix- 
ing the widths a priori, since the dynamic range for the various 
parameters is too high. Such optimization of the proposal distribu- 
tions before the main run speeds up the process, wh ile preserving 
the ergodicity of the MCMC algorithm tGelfand & Sahu. 1994) . We 
set T to 0.9, which results in points A x'^ ~ 1 away from the best 
solution to be accepted with 30% probability. Note however, that 
since increasing the temperature increases the acceptance rate as 
well (preserving the Ax^/r ratio) and vice versa, in principle the 
temperatiffe and g-width are degenerate, and therefore the exact 
value of the temperature chosen is not crucial as it will be compen- 
sated for in the above width adjustment. For the overall length of 
the chain, we need to find the length which both samples the poste- 
rior probability distribution well and at the same time (for practical 
reasons) is not much longer than what is just needed to accomplish 
this. We use 600,000 iterations, since our experience shows that the 
posterior probability around the best solution is well-sampled by 
this time. 

This procedure still leaves us with little sensitivity to highly 
disjointed solutions of equal goodness-of-fit, although this is not 
of concern here as such solutions will most likely converge onto 
unphyiscal values. 

Since we start at a high-probability part of the parameter 
space, the 'burn-in' period is less well defined than when starting 
at a random position, and we therefore do not formaly subtract a 
'bum-in' part of the chain. However, in deriving the best-fit param- 
eters and errors below we apply a cut of Ax^ = 5 from the minimum 
solution, which isolates the region of interest and thus has the 
same function as the 'burn-in' removal. Fig. lAll shows an example 
of the resulting chain in various projections. 

A2 Error estimates 

As discussed above, the MCMC procedure samples the likelihood 
surface, which is proportional to the posterior probability distribu- 
tion of interest. We use flat priors, where negative values are not 
accepted for any parameter In addition, upper limits are set for the 
temperature (100 K), optical depth (30), and 13 (3). The latter is nec- 
essary, since for many sources the SNR in the sub-mm data is very 
low, leading to the difference in overall x^s for unphysically large 
values of 13 to be small, and the chain can get stuck in such regions. 
Since 13 is believed, both on the oretical and observ ational grounds, 
to be in the range ^1-2 (see iDupac et aljr2003l and references 
therein), we believe our 0-3 prior is reasonable. The temperature 
limit is also necessary because of the above difference in SNR. The 
optical depth limit is only relevant for the few sources with higher 
r and faint near-lR/IRAC detections (e.g. Nl-040, Nl-064). These 
limits are imposed by returning the chain to the vicinity of its initial 
state in order to avoid getting stuck at the edge. We come back to 
the effects of the prior below. 

The chains obtained above allow for two distinct routes to ob- 
taining the probability distributions for each parameter The first is 
the easier straight marginalization of the given parameter over all 
the others. This is represented by: 

p{aj)daj = / p{ai)dai, (A3) 
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Figure Al. Here we demonstrate what the resulting chain looks like after 
the procedure described. The two left-most panels show examples of both 
uncorrelated and correlated parameters. A thinning factor of 30 has been 
applied. Throughout, we plot all points with ^ (x^ ■ +5). The top 
right panel shows a section of the chain demonstrating the evolution of the 
PAH amplitude parameter The bottom middle and right panels show the 
mai'ginalized distribution for the /3 and parameters (with linear y-axes). 
Note the non-Gaussian shapes. 



which in practice is just counting the number of times the chain 
visits a particular bin of values for the given parameter aj . How- 
ever, since we record the value for each chain link, we can also 
directly obtain the likelihood distribution (see ea-n lA2> for the pa- 
rameter. In p ractice, a simplified form of this is the mean likelihood 
distribution iLewis & B ridle 2002). The idea is to calculate the av- 
erage x^ for each bin and then compare this with the minimum 
achieved by the chain as 

p{aj) oc exp[-(xj - X?nin)], (A4) 

where x| is the mean x^ in the bin. 

In Fig. |A2l we show the distributions for all parameters ob- 
tained using both methods for the whole sample. Note that since 
we have applied a x^ cut (see above), secondary features in the 
marginalized distributions, such as blended peaks and tails, are still 
fairly likely. For example, in the case of Nl-101, Td ~ 25 K is the 
preferred solution, but Td ~ 40 K is still quite likely. The uneven 
error-bars indicate this. Note that in Table |2| for simplicity, we 
quote the average error only, but indicate the different error-bars 
in the figures. However, due to its definition, applying a x^ cut 
when estimating the mean likelihood tends to flatten the distribu- 
tions. Therefore in that case, the most 'peaked' distributions pos- 
sible are obtained if the whole of the chains are used. This means 
that parameters poorly constrained by the x^ are clearly visible in 
the likelihood curves (e.g. Lpah for the highest-^ sources Nl-064, 
andNl-078). 

The best-fit values we use are merely the peaks of the above 
distributions. To obtain the uncertainties on those we define a de- 
sired confidence level and, starting from the peak of the distribu- 
tion, move to incrementally lower equal probability bins on each 
side of it until the area covered is equal to the total area times the 
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Figure A2. Here we compare the marginalized (solid curve) and mean like- 
lihood (dot-dash curve) derived probability distributions for the parameters 
of interest. The distributions are normalized and the plots scaled to just fit 
the distributions as shown. The sources here and in the following figures 
are in order of increasing Lir. We use the marginalized distribution for our 
error estimates, where the best-fit value is given by the solid line (simply 
the peak of the distribution), and the 68% confidence limits are given by the 
dotted lines. 

desired confidence level. Thus unsymmetric errorbars are indica- 
tive of asymmetric probability distributions. 

In Fig.lA2l we see that in general the marginalized and mean 
likelihood distributions agree reasonably well with each other. 
However, there are two instances where they disagree substantially. 
One is for (3, where the likelihood tends toward unphysically high 
values. The origins of this are discussed in Appendix C. Our im- 
posed prior confines the marginalized distribution to more reason- 
able values. We also find significant differences between the two 
approaches for tv- Here we are again clearly affected by our choice 
of prior. Fig. IA3I shows how our uniform prior-based distribution 
transforms into the mean likelihood distribution by application of 
the Jeffreys prior (see e.g. Gregory et al. 2005). Note that the off- 
sets in the PAH luminosity and stellar mass values are also due 
to this difference. For the highest-z sources, Nl-064 and Nl-078, 
the marginalized distributions result in essentially unconstrained 
TV (no clear peak in the distributions), thus the values listed in 
Table 2 are fairly meaningless. For these two cases, we find the 
mean-likelihood distribution to be more indicative, in both cases a 
clear, although fairly broad, peak exists at ry ~ 5. We discuss the 
effects of this discrepancy as appropriate. 



APPENDIX B: PHOTOMETRIC REDSHIFTS 

The most common approach to photometric redshift estimates 
is template fitting. Its effectiveness is obviously strongly de- 
pendent on the templates used and the data available (see e.g. 
[Massarotti et al. 2001, for a discussion). 
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Figure A3. The effect of choice of prior on the tv probability distribution. 
Here we use the chain for Nl-041. 
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Figure Bl. The marginalized distributions for the photometric redshift esti- 
mates. The best-fit value is indicated by the sohd vertical line. Where spec- 
troscopic redshifts are available they are indicated by the dashed vertical 
line. Where spectroscopic redshifts are not available the best-fit photomet- 
ric redshift is indicated in brackets below the name of the source. 

Our multi-component SED model (described in Section IT2l 
effectively represents a large and flexible template library. There- 
fore, we here describe the results of fitting this model, while keep- 
ing the redshift a free parameter. The results we present use the full 
IR SED, however, we note that the near-IR and IRAC data are the 
most constraining in terms of redshift determination. This is par- 
ticularly true for z ^0.3 where the major PAH features leave the 
IRAC 8 fim band. With sufficiently good SNR, the approach also 
works at somewhat higher redshifts due to the 1.6 /im stellar peak 
being probed by the near-IR data (although in our case the IRAC 
data is not sufficiently sensitive to do this reliably at z ~ 1). We use 
model spectra which are generated on a logarithmic grid in wave- 
length, where A log A = 0.005. This leads to a minimum redshift 
resolution (Az/{1 + z)) of 0.01, which is adequate for our pur- 
poses. 

In order to adapt our earlier model-fitting to this problem, a few 
modifications are needed. First, our results in Section 4.4 sug- 
gest that /3 2 adequately described the far-IR/sub-mm emission 
of our galaxies. Therefore, here we fix /9 = 2. Next, to avoid some 
clearly unphysical solutions, we restrict the code to the range - 
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0.1 < log(^star/^PAH) < 0.9. This Condition is fully satisfied by 
all our sources when the spectroscopic redshifts are used (except 
Nl-040, where PAH emission is apparently completely absent), and 
moreover is required if the rest-frame sources a re to agree with our 
simulated IRAC colour-colour plot in lSaiina'et a l. ( 2005). Lastly, 
due to the T/(l + z) degeneracy, the temperature parameter we 
vary here is T / (1 + jz). This last step removes the strongest cor- 
relation in the resulting chains. However, we find that the results 
are subject to additional degeneracies of redshift with other param- 
eters, including: a negative correlation with the stellar amplitude; a 
positive correlation with the PAH amplitude; and a negative corre- 
lation with optical depth. These are clearly not independent. 
Fig. lBll shows the marginalized distributions we obtain for photo- 
metric redshift fits. Note that the mean likelihood distributions here 
are poorly constrained, most likely due to the above degeneracies. 
Overall, the agreement between our MCMC estimates and the spec- 
troscopic redshifts (where available) is quite good. The worst case 
is Nl-040 where the near-IR/IRAC data have very low SNR and 
moreover the source is apparently devoid of PAH emission, thus 
the constraint comes entirely from the far-IR/sub- mm which is not 
very a ccurate, as can be seen in this case (see also iAretxaea et alJ 
[200V)). We have left out Nl-064, whose high redshift (0.91) is in a 
different regime (the PAH features have completely left the IRAC 
bands) and therefore an accurate redshift cannot be obtained with 
the above method. Nl-078 is also believed to be at z ~ 1 (S03) and 
has an SED very reminiscent of Nl-064. Therefore we will assume 
it to be at 2 = 0.91 as well, although note that the uncertainty is 
~0.5. 

Fig. lBll shows that there is a bias in the model towards prefer- 
ring somewhat higher redshifts compared with the spectroscopic 
ones. The above degeneracies result in double-peaked distributions 
(e.g. Nl-OOl, Nl-007), and in some cases the 'wrong' peak is 
preferred (N 1-004, N 1-031). All of these issues suggest that the 
model, although working in general, is imperfect at this point. Im- 
provements are needed in minimizing the remaining degeneracies, 
which must involve including additional constraints. Two obvious 
steps are the use of the relation between absorption-derived and 
emission-derived dust mass, as discussed in Section 5.2 (this re- 
quires the effective sizes of galaxies), and the use of distance dim- 
ming. The latter for example effectively excludes the higher-z solu- 
tion for N 1-004, but not for N 1-031. For Nl-040, the distance dim- 
ming argument would in fact support the wrong solution, which 
was brought about largely by the imposition of 'typical' stars-to- 
PAH ratios. Therefore, the range in galaxy SED types needs to be 
known better and accounted for here before these steps can be im- 
plemented. 

After accounting for all these issues, we find that the average 
uncertainty of our estimated redshifts is (A[2/(l + z)\) =0.06. 

Thus in summary we adopt the following redshifts: Nl- 
013(0.46), Nl-045(0.18), Nl-068(0.16), Nl-077(0.20), and Nl- 
083(0. 31). Based on optical template fitting lRowan-Robinson et ^ 
<2005h derived the following redshifts: Nl-013(0.58), Nl- 
045(0.09), Nl-068(0.23), Nl-077(0.19), and Nl-083(0.51). As- 
suming these are correct, our results on average are uncertain to 
{A[z/(1 + 2:)])=0.07, which is slightly worse than the above 
quoted uncertainty as expected when comparing against other pho- 
tometric redshift estimates. For simplicity, we are going to ignore 
the uncertainty on the estimates and treat these redshifts as fixed. 



APPENDIX C: THE T - /3 RELATION 

A common approach, especially with a limited number of far- 
IR/sub-mm data points, is to fit a greybody function with charac- 
teristic temperature and emissivity. There are three main questions 
arising from the practice and interpretation of these fits: 

1) what degeneracies are present in the usual formulation, and 
therefore what are the optimal parameters to fit? 

2) given that this is only an approximation to the true shape of the 
SED, what effect does the spectral sampling (such as due to red- 
shift) have on the results obtained? 

3) how does one disentangle physical correlations from said degen- 
eracies? 

For further discussion, see alsolBlain et alJ (20031) and the al- 
ternative parameterization of Baugh et al. (2005). In essence, the 
function one is fitting has the form: = {v/uo)^ x , where 
By is the Planck function: 



B, = 



2h 



c2 exp{hu/kT) 



(CI) 



The parameters to determine here are Td , /3, and the nuisance 
parameter (or some equivalent normalization). We wish to dis- 
entangle any correlations between the three. To begin with, to un- 
derstand how degeneracies arise between the various parameters, 
we need to consider what it is that the code is actually fitting. Intu- 
itively, to first order this is the overall amplitude of the dust emis- 
sion, the position of the peak, and the slope in the Rayleigh- Jeans 
tail. 

The first of these functional parameters is the bolometric dust 
flux, Fd, which is given by 



= / .Udu = 



2h 



f-V 



r(4 + /3)C(4 + /3), 



(C2) 



where F is the complete gamma function, and ( is the Riemann 
zeta function*. Rearranging the above and substituting for l/i^o 
in fi, removes the bulk of the degeneracies of the normalization 
parameter (which is now taken to be _Fd). 

Nevertheless, the temperature and /3 are still correlated. As 
stated above, the fitting process is also concerned with the position 
of the peak in the emission. By solving d/^ /di^ = we obtain: 



r, _ ( fc'^pcak \ 



(C3) 



This equation provides a good fit to the observed correlation at 
higher values of (3 (see Fig. CI), but fails as 13 approaches 1. The 
reason for this is obvious - apart from the peak, we also need to 
match the Rayleigh Jeans tail. The sub-mm data do not allow [3 to 
become arbitrarily shallow, although it can be compensated by in- 
creasing the temperature (which pulls the Planck function the other 
way). The constraint coming from the Rayleigh-Jeans tail can be 
parametrized with the sub-mm spectral index, Qsubmm which is re- 
lated to [3, but in the Rayleigh-Jeans approximation the temperature 
dependence falls out leaving: 



21og(850/450) 



(C4) 



Looking at our data (when the SNR at 450 /xm was high), we find 
typical values of asubmm ~ 5, which means g~ 1.5 from ea-n.lC4l 



* Note that (J(4 + /3) is ~ 1 for any positive value of (3, so does not affect 
the result much. The quantity logr(4 + /3) is closely approximated by 



0.69(1+0.69/3), which we use. 
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discussed above) as the peak of the SED shifts, the combined /3, T 
distribution moves as well. But in the Dupac et al. sample the short- 
est wavelength is 100 /^m, thus whenever the peak is shortward of 
that (corresponding to ~ 30 K for l3 ~ 2), the sub-mm index be- 
gins to dominates the fit, flattening the distribution. At the other ex- 
treme for example, as the peak reaches the coldest observed values 
(corresponding to peaks at ~ 200 /im), we are sampling a distribu- 
tion where r'^15Kfor/3~2 (see Fig. CI). Therefore, it appears 
that this correlation is easily explained as a sequence of progres- 
sively colder/hotter environments. The apparent l3-T correlation 
is merely a combination of the intrinsic correlation between Td and 
/3 for individual sources and limited wavelength coverage. That be- 
ing said, a physical correlation here is to be expected from the dif- 
ferent optical properties of grains at different temperatures. But it 
is unclear that the single greybody approach is capable of testing 
this adequately. 

This paper has been typeset from a TjiX/ KTgX file prepared by the 
author. 



Figure CI. Here we compare the relation in ea-n lC5l and the observed de- 
generacy between the /3 and temperature parameters. The points are from 
the fit to Nl-024. The dashed lines show the relation in eq-n IC3l where 
the location of the peak moves in increments of 20 fim from 100 /^m to 
200 ^m (top to bottom). Note the good agreement at higher /3's and the de- 
viation at lower-valued /3's. The dotted line is the expectation of eQ-n lC4l 
with Cfsubmm=5 (consistent with observations). The solid lines are the 
combined constraint (eg-n lCSl for peak locations of 160 /^m and 180/tm 
respectively - the observed distribution is now reproduced much better. If 
not treated carefully, this parameter correlation, which exists in the fits for 
each source, can be misinterpreted as a physical coiTelation between the 
derived (/3,T(j) pairs for a sample of sources. 



which is in fact often adopted for sub-mm sources. The joint con- 
straint on the spectral index leads to: 

/850\ ^ 

+ asubmm-ilOg ] -D.(C5) 

In Fig. CI, we show how the above two constraints act together and 
therefore how the available data drive the f3 and temperature values 
found. For example, in data without well determined peak position, 
assuming a low value of /3 (say, 1-1.5) will automatically lead to 
the conclusion of hotter temperatures. Conversely, in data with cool 
temperatures and well determined peak, but poorly sampled sub- 
mm data, the conclusion will always be that /? is high (even > 2 if 
the fit is allowed to go there). Better sampling (which could mean 
better SNR or more spectral points) in the sub-mm (relative to the 
peak) would likely change that un-physical conclusion. Conversely, 
as the redshift is increased, our data increasingly sample the peak 
rather than the tail of the thermal emission, and this has the same 
effect. 

This degeneracy is empirically well known tel ain et alj2003h . 
Apart from this degeneracy arising from the functional form used, 
there has also been suggestions that a physical inverse relationship 
exists as well. In particular, this was seen when the best-fit tem- 
perature and p values for different Milky Way environments were 
compared by Duoac et al. ( 2003). The relative homogeneity of our 
sample (see Fig. l4.4t does not allow us to firmly support or reject 
this relationship (although the apparently warmest source, Nl-078, 
supports it). A wider range of galaxy types with well sampled SEDs 
would be needed to address this firmly. However, we note that (as 
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